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Abstract 



We analyze the Optimal Channel Network model for river networks using 
both analytical and numerical approaches. This is a lattice model in which a 
functional describing the dissipated energy is introduced and minimized in or- 
der to find the optimal configurations. The fractal character of river networks 
is reflected in the power law behaviour of various quantities characterising the 
morphology of the basin. In the context of a finite size scaling Ansatz, the 
exponents describing the power law behaviour are calculated exactly and show 
mean field behaviour, except for two limiting values of a parameter characteriz- 
ing the dissipated energy, for which the system belongs to different universality 
classes. Two modified versions of the model, incorporating quenched disor- 
der are considered: the first simulates heterogeneities in the local properties 
of the soil, the second considers the effects of a non-uniform rainfall. In the 
region of mean field behaviour, the model is shown to be robust to both kinds 
of perturbations. In the two limiting cases the random rainfall is still irrele- 
vant, whereas the heterogeneity in the soil properties leads to new universality 
classes. Results of a numerical analysis of the model are reported that confirm 
and complement the theoretical analysis of the global minimum. The statistics 
of the local minima are found to more strongly resemble observational data on 
real rivers. 
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1 Introduction 



Experimental observations on river networks have shown clear evidence of their frac- 
tal character. Data from many basins with different geological features have been 
analyzed, and have shown power law behaviour of the probability distributions of 
several quantities describing the morphology of the river basin [|l], H, H, H over a wide 
range of scales. 

Several statistical models have been proposed ^ ||, but a complete the- 

oretical understanding is still lacking. Recently a simple lattice model derived from 



an energy minimization principle has been proposed [TT|, [12], [T^ that in spite of its 
simplicity, seems to reproduce many features of natural river networks. 



Numerical investigations of the model have been performed p, |T^, 0, |1^ with 



different geometrical constraints on the form of the basin. Furthermore, the model has 



been analyzed within the framework of a finite size scaling Ansatz [^] . In the present 
paper the so called Optimal Channel Network model [|ll|, |l^ is studied analytically 
and exact results are obtained. 

In addition, generalized models, taking into account the presence of quenched disorder 
are considered. Randomness is introduced in two different ways: one modelling the 
inhomogeneity of the soil and the other non-uniformity in the rainfall. Analytical 
results are given in these two cases, too. 

Results of numerical simulation in the last part of the paper, confirm and complement 
the analytic predictions. 

In Section 2 we describe the lattice model and derive the scaling laws. The rela- 
tionship between exponents are also derived. The exponents characterizing the power 
law distributions of drained areas and mainstream lengths are expressed in terms of a 
single independent exponent: the wandering exponent in the self-affine case and the 
fractal dimension for the self-similar basin. The section ends with the definition of 
the Optimal Channel Network model and with a short discussion of the underlying 
minimization principle. 

Section 3 is entirely devoted to a analytical study of the model on a simple fractal, 
the Peano basin. The solution is given exactly and is used to give bounds in the 
subsequent Section. The distributions of areas and lengths are evaluated and shown 
to exhibit the form predicted by the finite size scaling Ansatz. 

In Section 4, analytical results in the absence of disorder are derived. The model is 
shown to exhibit three distinct universality classes for different values of a parameter 
characterizing the dissipated energy. 

Heterogeneities in the soil properties and random rainfall are considered in the gen- 
eralized models of Section 5. Analytical results for these cases are also deduced. 
Numerical results pertaining to the search for the global minimum of the dissipated 
energy with a simulated annealing algorithm are given in Section 6 and numerical 
results for the statistics of the local minima are given in section 7. 
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2 Definitions and Derivation of Scaling Laws 

A river basin is described by a scalar field of elevations. Drainage directions are 
identified by steepest descent, i.e. by the largest local decrease of the elevation field. 
The presence of lakes has not been considered, i.e. from each point the water can 
fiow to one of the nearest neighbours. This corresponds to have all lakes saturated. 
Within this context, a river network can be represented by an oriented spanning-tree 
over a two dimensional lattice of arbitrary size and shape, in which oriented links 
(one coming out from each site of the lattice) correspond to drainage directions. 
We will consider spanning trees rooted in a corner of a L x L square lattice (the 
root will correspond to the outlet). Site i is upstream respect to site j if there is an 
oriented path from i to j. To each site i of the lattice, we associate a local injection 
of mass Ti (the average annual rainfall in the site i). The fiow Ai, also referred to as 
the accumulated area can thus be defined as the sum of the injections over all the 
points upstream of site i (site i included). 
The variables are thus related by 

Ai = J2'^i,j^j + ri, (1) 
j 

where Wij is 1 if site j is upstream with respect to site i and is a nearest neighbour 
of it and otherwise. The local injection rj is commonly assumed to be homogeneous 
and identically equal to 1. 

In natural basins these areas can be investigated through an experimental analysis of 
digital elevation maps (DEM's) 0. See figure la for an example. 

The upstream length relative to a site is defined as the length of the stream 
obtained starting from the site and repeatedly moving in the upstream direction 
towards the nearest-neighbour with biggest area A (the one leading to the outlet is 
excluded, since it is a downstream site), until a source is reached i.e. a site with 
no incoming links (see Fig. lb). If two or more equal areas are encountered, one is 
randomly selected. 

For a given tree, one may consider the probability distribution of the following quan- 
tities: for a lattice of given linear size L we will call p{a, L) the probability density 
of accumulated areas a and 7r(/, L) the probability density of the upstream lengths /. 
These represent the fraction of sites with area a or stream length / respectively. We 
will also consider the integrated probability distributions P{a,L), the probability to 
find an accumulated area bigger then a and n(/, L), i.e. the probability to have a site 
with an upstream length bigger than I. 

Both these probability distributions, here defined in the simple case of the lattice 
model, were originally introduced to describe real rivers and experimentally found to 
scale as power laws leading to the formulation of a finite size scaling ansatz[p!7[] 

p{a,L)=a-^f(—), (2) 
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7r(/,L) = r'^^(l), (3) 

where f{x) and g{x) are scaling functions accounting for finite size effects and ac and 
Ic are the characteristic area and length respectively. The functions f{x) and g{x) are 
postulated to have the following properties: when x ^ oo they go to zero sufficiently 
fast to ensure normalization; when x they tend to a constant, to yield simple 
power law behaviour of the probability distributions in the large size limit. This also 
implies that r and ip are bigger than one. 
The characteristic area and length are postulated to scale as 

ac ~ , (4) 

Ic ~ L''' . (5) 

In river basins, anisotropies are always present due to a non-zero average slope of 
the landscape and the presence of gravity. Thus, one has to distinguish between a 
typical longitudinal length L and a typical perpendicular one L± (these two length 
are measured along the two principal axis of inertia), which scale as 

L± = L'', (6) 

giving ac ~ L^~^^ , i.e. (f = 1 + H. H is known as the Hurst exponent and of course 
< if < 1. In what follow, for the sake of simplicity, we consider basins of square 
shape; the above relations then, refers to dimension of a generic sub-basin inside the 
bigger one, whose dimensions are fixed from outside. 

The exponent thus corresponds to ip = 1 + H . The di exponent, characterizing 
the typical length, can be assumed to be the fractal dimension of a stream (for fractal 
river networks, each rivulet going from any site to the outlet is a fractal with the 
same fractal dimension) , and is such that 1 < di < 1 + H . The bounds correspond 
to a straight line and a space-filling behaviour. For self-affine river basins we expect 
di = 1 and H < 1, whereas, when H = 1 then di > 1 in the self-similar case. 

For the same quantities, the integrated probability distributions can be analo- 
gously written: 

P(a,L) = a^-F(^), (7) 



U{l,L) = l^-^G 
which follow from (13) and (|3D with 



I 



Fix) = X--' / dy y-^fiy) , (9) 
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G{x) = x^-' / dy y-^g{y) , (10) 

where sums over variable y have been replaced by integrals in large L limit. 
Prom the above definitions and the properties of / it easily follows that 

(a") ~ L^^+H){n-r+l) ^ 

for any n > r — 1, while (a") ~ const, if n < r — 1. Note that both a and I have a 
lower cutoff which is one. 

This equation, evaluated for n = 1 gives for the average area, 

(a) ~ L(i+^)(2-r) _ (12) 

The mean accumulated area (a) can be easily shown to be equal to the distance from a 
site to the outlet, averaged over all sites. In effect, in the sum over all the downstreams 
(the rivulets going from each site to the outlet), the number of times each bond (of 
unit length) appears exactly equals the accumulated area of the associated site. Thus 
summing over all Ai is equivalent to a sum over all the downstream lengths: 

(q) — {I downstream) ■ {^^) 

{^downstream) Can be evaluated replacing in the sum the distance of each point from 
the outlet measured along the stream with the corresponding euclidean distance d{x) 
to the power df. 

{^downstream} "To' ^ ^ Idownstreamj-^) "To' ^ ! d(^x) ' ~ Z/ ' (^4) 
^ X ^ X 

This fact is general and the argument we used does not need the knowledge of the 
downstream length distribution. This distribution however can be explicitly derived 
at least in the case of directed networks. We call directed those networks such that 
each oriented link has positive projection along the diagonal oriented towards the 
outlet. 

The reason to introduce this class of networks is that river basins often have a quasi- 
directed character, due to the fact that they are typically grown on a slope which 
gives a preferred direction to the flow. Moreover directed trees are much more simple 
to handle analytically than "undirected" ones. 

For such trees consider the set of 2L diagonals orthogonal to the one passing through 
the outlet: the downstream length is the same for all the points on the same diagonal. 
Thus the number of points at a given distance to the outlet can be easily seen to be: 

number of point at distance / = < ^ , \ _\\^^ or (^^) 

thus the probability density for the downstream lengths has the form of a power law 
with exponent —1 times a scaling function of the argument IjL: 
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downstr earaij' I I f downstream^ (-^6) 



with 

fdownstreamix) = min(x^ 2x - X^) <X <2 (17) 

The first moment of this distribution again gives eq. ([I^ with di = 1 which is the 
expected result for the fractal dimension of a directed tree. This resuh, together with 
|14| suggest that downstream length distribution might have the scaling form: 

downstr eamij' I / fdownstreami j j ) (-^^) 



for the general case. 

Equations (|13]) and (|l^ lead to the following expression for the average area: 

(a) ~ L^' . (19) 



From eq.(19) and eq.(12) we get the scaling relation 



l + ^-p^. (20) 

A well known hydrological law, Hack's law |T^, relates the length of the longest 
stream / in the drained area to the drained area a of the basin: 

/~a^ (21) 

The accepted value of /i is /i = 0.57 ± 0.06 |T^, whose difference from the 

Euclidean value 0.5 lead to the first suggestion of the fractal nature of rivers ||l[. 
From equations (^) and (^ it follows that 

h = . (22) 

1 + H ^ ^ 

Together with vr and p one can define the conditional probability vf (/ | a) of finding 
a main stream with length / in a basin with accumulated area a. The simplest scenario 



is that (^) still holds and 7r(/ | a) is a sharply peaked function of / with respect 



to a, i.e. there exists a well-defined constraint between lengths and areas, 

7f (/ I a) = 6{l - a^) , (23) 



or more generally 



nil \ a) = 1-^1^]. (24) 



For the probability density vr, p and vr the following consistency equation must hold: 
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7r(/,L) = / da n{l \ a)p{a, L) (25) 



that gives, in the large L hmit 

{ij-l)di = {r-l){l + H) (26) 

relating the exponents in the distribution of lengths and in the distribution of accu- 
mulated areas. 

The scaling relations ( ^0]) and (pGf) can be expressed in a simpler form, observing 
that both r and ijj depend on di and H only in the combination j^^^ = h, where h 
is the parameter appearing in Hack's law (pi]). Thus 



r = 2-h, (27) 
* = i. (28) 

The exponents characterizing the distributions of accumulated areas and upstream 
length are thus related by the simple expression: 

r = 2 - i (29) 

For self-affine river basins 

H <1 ,di = l, (30) 
and all exponents can be expressed in terms of the Hurst exponent H, giving: 



1 + 2H 

= l + H- (32) 

while in the self-similar case 



yielding: 



H = 1 ,di>l, (33) 



r = 2-|, (34) 
^ = J- (35) 
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Note that in both cases, r < 3/2. The equahty holds only when H 
corresponds to the mean field situation |2^. Likewise, h >l/2. 



di = 1 which 



A recently formulated lattice model [|Tl|, [12], ^ based on a minimization principle, 
seems to reproduce quite well the main characteristics of river networks. In this model, 
a rule for selecting particular configurations in the space of spanning trees is given. 
The "right" configurations, called Optimal Channel Networks (OCN) are obtained on 
minimizing a dissipated energy, written as: 



E = Y. hAz{t)Q, 



(36) 



where Qi is the fiow rate (the mean annual discharge) in the bond coming out from 
the site i, Az{i) is the height drop along the drainage direction, and ki characterizes 
the local soil properties such as the erodability. It will be taken to be equal to one 
for each bond for homogeneous river networks. 

Given a field of elevations, drainage directions are usually identified by steepest de- 
scent, i.e. by the largest downward gradient Vz{i) of the scalar field z{i). This allows 
us to obtain another expression for the dissipated energy on adding the following 
considerations: 
i) in the case of uniform rainfall in time and space 

Qi ~ 

a) experimental observations on rivers suggest an empirical relationship between the 
fiow rate and the drop in elevation: 

Az{z) ~ Q]-' 
with a numerical value around 0.5 for 7. 

Thus one obtains, apart from a multiplicative constant, the alternative expression of 
equation (|36|) : 



(37) 



which was first proposed by Rinaldo et al. |Tl], |T2|, |2^, and will be analyzed further 
in this paper. 



3 The Peano Basin 

The Peano basin is a deterministic space filling fractal on which exact calculations 
can be carried out 
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It has a spanning tree like structure not too dissimilar to 
that of real basins. The scaling laws for such a basin can be obtained exactly and 
the dissipated energy (|37|) can be estimated. The latter provides an upper bound for 
the minimum energy dissipated by an OCN and will be a crucial ingredient for the 
calculations to follow. 
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The Peano basin is obtained as follows: start with an oriented link; at step 2 such 
a link generates four new links, two resulting from the subdivision in half of the 
old link and preserving its orientation, and the other two having a common extreme 
in the middle point of the old link and both oriented toward it (see Fig. 2a). At 
each successive step, for each link four new oriented links are substituted in the way 
previously described. After T steps the fractal has Nj- = points (excluding the 
outlet) and it can be mapped on a square lattice of size L = 2^ with bonds connecting 
first and second neighbours to form a spanning tree. 

We can associate with each site i of the Peano basin (iterated until step T) an 
area Ai{T) as in the previously defined lattice model of river networks. Let Vt 
denote the set of distinct values assumed by the variable Ai at step T. It can be 
easily checked that Vt contains Vt-i and 2^~^ new distinct values, appearing for the 
first time. Thus Vt contains distinct numbers. Let us call A = Ur=o and a„ the 
increasing sequence of numbers in A (the distinct values of Ai that can be generated 
iterating the construction). For such a sequence, the following rule holds 

an = S{Y,Ck{n)A') + l n = 0,l,... (38) 

k 

where the Cfc(n) are the coefficients of the binary expansion of n (n = I]fcCfc(n)2^). 
In the construction described in Fig. 2b, denote by Mj the number of sites i with 
Ai = ttn at step T. The following recursion relation then holds: 

Mj = 4-Mj-i-l T>t(a„) 

Mj = 1 T = t(a„) (39) 

Mj = T < t(a„) 

where t(an) is the first step in which an area with value a„ appears, and is given by 

. X _ J n = , . 

H«nj-| i + [iog2(n)] n>Q ^^^> 

where [■] is the integer part. 
Solving ( p9D one gets 

_ f , T < t{ar. ^ 



Mn - { 2^T-t{a„) _^ 1 ^ rjn ^ ^( J" \ (41) 



3 '3 

and thus all the a„ "born" at the same time step have the same probability Pridn) = 
p{a^,L = 2^) = Ml/NT 

Then the integrated distribution of areas P{Ai >an,L = 2^) assumes a very simple 
expression for a„ of the form 4* (one can easily check from (|38D that a2'-i = 4*), and 
is given by 



P{A,>a = 4:',L = 2')=a'-^F( ) (42) 



10 



having the form (|^) with 

r = 3/2 , H = 1 (43) 

and 

F{x) = ^{l~x) 0<a:<l. (44) 

and F{x) = when x > 1. 

Equation ( ^2]) is obtained on observing that P{Ai > a = 4* , L = 2^) = Z]n=2* PricLn) 
depends on n only through t(a„), allowing one to replace the sum over n with a sum 
over the steps s. Moreover, for each step s > there are 2^*"^ areas with the same 
t{0'n) = -s. Thus 



P{A>a = 4\L = 2^)= !4(^-^) + I " ^ = 

/ ^ \ ^ (45) 

- \2-*fl-22(*-^)) = ia-Ul-^l. 



Similarly, choosing / of the form / = 2* and observing that at step T the sites with 
upstream length greater than or equal to 2* are the ones in which the accumulated 
area exceeds 4*, we find that 

n(Z >2\L = 2^) = 1^-^g(^-^^ (46) 
which is of the form (^ with 

^ = 2, di = l (47) 

and 

G{x) = \{l-x'). (48) 

In the following Section we will need some estimates of the mean value of for 
a Peano basin of size L = 2*. 

-1 oo 

(^') = 7iE^7 = Ep-KK"- (49) 

^ i n=0 

From the expression for a„, it follows that 

^ ■ 4*("") < a„ < 4*('^") (50) 
(t(a„) is the "birth-time" for a„), giving 

^)'a(L = 2^)< (A^) <«(L = 2^), (51) 
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where 



„(i.2n.|,-L^|(24(.-..,l)!_i4... ,52) 
Performing the summation one gets in the large size hmit 

r 1(1 + ^32^) , 7<i/2 

«~ 31^2 log^ ' 7 = 1/2 (53) 
[ T^L^^'^ , 7 > 1/2 

From equation (|5TD, one gets 

{const 7 < 1/2 
logL 7 = 1/2 (54) 
7 > 1/2 

which will be essential to obtain an energy bound for the lattice OCN model. 

The scaling exponents for the Peano basin can also be obtained by a renormaliza- 
tion group argument. Let us consider, for example, the scaling of the accumulated 
areas. 

The self-similar structure of the Peano basin suggests a natural "decimation" proce- 
dure. The idea is the following: consider the equations relating areas at time step 
T; then, eliminate the variables related to the sites that are not present at time step 
T — 1. This leads to an effective equation describing the same physics on a tree scaled 
down by a factor 2. 

For the sake of simplicity let us consider the Peano basin at the second step of 
iteration. In Fig. 3 let A^^^ denote the variables related to sites that are present at 
step T = 1 and B'^'> denote the ones that will be eliminated by decimation. The 
upper label refers to the step. In what follows it will be useful to write the equations 
in terms of A^^^ = — 1 and = — 1. The areas at step T = 2 are related 
by: 

Af^ = 3 • ^ + 3 , 

fif ) = ii^) + 2 ■ 5f + 3 , (55) 

^0 

Elimination of the B^^ leads to 



= . 



Af^ = 3 ■ + 12 . (56) 
At time step T = 1 the relation between areas is straightforward: 

iji) = 3 ■ + 3 . (57) 
Equations (^Bj) and ( |^ are the same if 
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i(T+i)=4if), (58) 

i.e. 



(Af - 1) = 4 (Af ) - 1) . (59) 

Denoting by n^'^'^^^a) the number of sites with area greater then a at step T + 1, one 
can easily observe that the number of decimated sites with A > a is half of the total 
number of sites with A > a 



n^^^'\a) = 2-n^^\a/4:). (60) 

Noting that the total number of sites at step T, iVy = 4"^, the integrated proba- 
bihty P(Af +1) >a) = ^ is given by 

P(Af+i)>a) = 6P(Af)>a/4) (61) 

with 

n(^+^(a)/4^+i _ 2/4^+1 _ 1 
~ n(^)(a)/4^ ~ ~T/4^ ~ 2 " ^^^^ 

The general solution of equation (|6T|) is of the form 

P{Al > a) (X a^-^ (63) 

apart from a superimposed oscillating term given by a periodic function of log(a) with 

period 1/4 pSf. 

From equations ( |6l] , |6^ ,|6^) 

r = l. (64) 

The same argument can be repeated for the distribution of stream lengths, recovering 
the ip exponent. 



4 Analytical Results: homogeneous case 

We now proceed to an analysis of the characteristics of the global minimum of the 
functional E for 7 in the range [0, 1] for the homogeneous case. 

Let us consider first the two limiting cases 7 = and 7 = 1. If we call /j 
the weighted length of the stream connecting the i-th site to the outlet (calculated 
assigning to each bond a length ki), it is straightforward to show that 

Y,kiA = Y,h. (65) 

i i 
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In effect, denoting with DS{i) {US{i)) the set of points downstream (upstream) with 
respect to the point i and observing that A{i) equals the number of points in the set 

US{i) one gets: h = T,i EjeD5(i) = T,j Y^jausii) kj = T,i hAi . 
The minimization of the energy dissipation for 7 = 1 thus corresponds to the mini- 
mization of the weighted path connecting every site to the outlet independent of each 
other. 

The 7 = case, on the other hand, corresponds to the minimization of the total 
weighted length of the tree: 



E = Y.h. (66) 

i 

In the homogeneous case fcj = 1, Vz, which leads to a high degeneracy for both 7 = 
and 1. 

Indeed for 7 = 0, each configuration has the same energy (each spanning tree on a 
L X L square lattice has — 1 links). For 7 = 1, the minimum of the energy is 
realized on a large subclass of spanning trees, namely all the directed trees, in which 
each link has a positive projection along the diagonal oriented towards the outlet. 
For the values of 7 G (0, 1) there is a competition between both mechanisms breaking 
the degeneracy and making the search for the global minimum a less trivial problem. 

The 7 = 1 case gives a minimum energy ~ L^. This can be derived observing 
that all points on a diagonal orthogonal to the one passing through the outlet have 
the same distance from the outlet. Then: E = Ylk=i ^(^ + 1) + Sfc=L^ k{2L — 1 — k) = 
L^(L — 1) ~ L^. Thus the value of the energy functional is the same for each directed 
network and corresponds to the Scheiddeger model of river networks - all directed 
trees are equally probable. Such a model can be mapped into a model of mass 
aggregation with injection that has been exactly solved by Takayasu et al. [^, ^ 
The corresponding exponents are: 



r = 4/3 , 7/; = 3/2 , if = 1/2 , rf, = 1 , /i = 2/3 . (67) 

These exponents follow easily from the result H = 1/2 and from our scaling solution 
of sec. 2. The result H = 1/2 can be deduced with a simple argument: since all 
directed trees are equally probable, each stream behaves like a single random walk 
in the direction perpendicular to the diagonal through the outlet. This implies that 
its perpendicular wandering is L±_ ~ L^/^. Comparing with equation (^ one gets 
H = 1/2 . 

The 7 = case leads to the same energy E L"^ for each network, thus reducing to 
the problem of random two dimensional spanning trees, whose geometrical properties 
have been calculated on a square lattice The results, in our notation, are: 



T = 11/8 , 7/- = 8/5 , if = 1 , = 5/4 , /i = 5/8 . (68) 
Both (^) and (|68D are consistent with the scaling relations (^), (|^) and (^9]) . 
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We now extend our analysis to the whole range 7 G [0, 1]. We will rigorously show 
30| that, in the thermodynamic limit, the global minimum in the space S of all the 



spanning trees of the functional E{^,T) = Y^iAiiT)^ scales as 

min E(7, T) ~ max(L^ (69) 

for all 7 G [0, 1] . 

Since -E(7, T) is an increasing function of 7 and it is equal to for 7 = 0, for 7 > 
it is obvious that 

E{-i,T)>L\ (70) 
We now observe that the sum over all the sites can be performed in two steps: 

2L-1 

E{ia)= E T.MTV (71) 

n=l 

where P„ are the lines orthogonal to the diagonal passing through the outlet, that 
we will enumerate starting from the corner farthest from the outlet (see Fig. 4). For 
directed spanning trees one can observe that the sum of the areas in a given line 
is independent of the particular tree, and is: 



k{k + l)/2 k<L 

- Sd{2L - 1 - k) L + l<k<2L-l ^ ' 



where 5^(0) = 0. This quantity only increases on generalizing to generic "undirected" 
trees and thus for an arbitrary spanning tree, 

S{k,T)= M'^) > Sd{k) . (73) 
Let us observe that for = 0, (L — 1) : 

S{k, T) + S{2L -l-k,r)> Sd{k) + Sd{2L -l~k)=L\ (74) 

making it convenient to perform the summation in (^) over "pairs" of lines. To get 
a lower bound for E we need a further inequality: for every set F 

that follows easily from Schwartz inequality, since Ai > 1 and < 7 < 1. Now, using 



T|), (0) and (ffSD we can write 



2L-1 L-l L-l 

Eh,T) = T.T.M'ry = i: E ^7>E( E Mr))' 

n=l i^Vn n=0 i eCDnUDJ n=0 i eCD„ui3„) 
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= [Sin, T) + S{2L - 1 - T)\' > L^^ = ^ (75) 

n=0 ra=0 



where Vn = V{2L-i-n)- The equahty in the last inequahty holds for (izrec tec? networks. 
We can thus write 

E{-f,r) > 1^+^-^ . (77) 
Equation (|77D together with the (^O]) gives the lower bound 

E(7,T) >max(L2,Li+27)^ ^yg) 

that holds for every T E S and thus also for the minimum over T. Using the results 
of the previous Section, we can exhibit a tree on which the bound is realized. In effect, 
the Peano network can be mapped on a square lattice only considering links between 
first and second nearest neighbours, but (|78|) can be analogously obtained for such 
a lattice on rearranging in an opportune way the summation in equation ([711). If we 
call Tp the spanning tree given by the Peano basin we know from equation (|5^) that, 
except for logarithmic corrections for 7 = 1/2, E{'y,Tp) ~ max(L^, L^'^^'^). Thus 



mm£;(7,T) <E(7,Tp), (79) 

and we get equation (|5UD. 

We now proceed to the calculation of the scaling exponents. 
For a directed path, from equation (pB]), 

(a) ~ L . (80) 
For generical "undirected" networks, let us write as in ([ToD 

(a) ~ L'^' (81) 

{di could be be somewhat bigger than one if one assumes a "quasi-directed" be- 
haviour). 

Using eq. ([TT|) with n = 7 and the above result on the scaling of energy: E = 
L'^{a') ~ L^+^T, one gets 



27-1 = (l + i7)(7-r + l) (82) 

holding for 7 > r — 1 . 

Equations ( p2D together with the scaling relation (^) can be solved with respect to 
r and if, and gives, for 7 > 1/2 

3(l-7)+(rfi-l)(l+7) 

^1-7)+^.-! (83) 
1-7 
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(the constraint 7 > r — 1 become 7 > 1/2, independently of di). Thus, if H < 1, for 
any •y < 1 di = 1, yielding 



4) 
of 



1, in which case (|8^) does not hold. 



r = 3/2 , H = 1 , 4j = 2 , di = 1 , h = l/2 

for 7 G (1/2, 1). The exponents are the same as in the mean field theory ||26 
the Scheidegger model and the same as in the Peano case |^ 
Note that equations ([83| ) are meaningless for 7 
consistent with the Scheiddeger results (|67D . 
When < 7 < 1/2, all we can say is that r > 1 + 7. However if = 1 for any 
7 G (0, 1) (i.e. the Optimal Channel is quasi — directed) then one might expect that 
H = 1 for all values of 7 and thus eq. ( p^ holds in the whole range 7 G (0, 1). This 
prediction is confirmed by results of numerical simulations. 



5 Analytical Results: heterogeneous case 

In this Section we focus our attention on the case in which some sort of quenched 
disorder is present in the basin. Two cases have been analyzed: 

i) random bonds, modelling heterogeneity in the local properties of the soil; 

a) random injection, modelling non-uniformity in the rainfall. 

In the first case, we will show that the energy can be bounded from above with 
the corresponding one in the absence of disorder. This will give in the large L limit 
an upper bound for the r exponent that we will see is realized for all the 7 G (1/2, 1). 
In the case of random rainfall, we will show that this kind of randomness does not 
affect the scaling behavior of the dissipated energy in the large L limit. All the 
analytical results found in Section 4 for the homogeneous case, being based on the 
energy estimate in the thermodynamic limit, can thus be extended to that case, giving 
the same values for the exponents. 

In the case of random bonds, we associate with each bond of the L x L square 
lattice (2L(L — 1) bonds) a quenched random variable k^, arbitrarily distributed such 
that (kh) = 1. The label b ranges over all the bonds of the lattice. The 2L{L — 1) 
variables are chosen independent of each other and identically distributed. 

In what follows b{i,T) will denote the label associated with the bond coming 
out from the site i on the tree T. Let T*{'~f) denote one of the trees on which the 
minimum of the energy E{'~f, T) is realized in the homogeneous case for a given value 
of 7, and S the set of all the spanning trees: 

E(7) = minE(7,r) = Y.MrV = Y.MT*{l)r . (85) 

i i 

Denoting by £"15(7) the minimum of the energy in the heterogeneous case, averaged 
over the quenched disorder, the following inequality holds: 
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En{^) = (min^A:,(,r)A(r)^)<(5:%,rn7))^.(^*(7))") 

i i 

= {T.hi^,r^M))Mr*{l)r = T.Mr*{i)r = e{^) . (86) 

i i 

The energy in the presence of this kind of disorder is thus bounded from above by 
the energy in the absence of disorder for any value of the 7 parameter. In the large 
size limit this result gives bounds on the scaling exponents. Equation ( |TI] ) evaluated 
for n = 7 gives 

which holds for any 7 > T£, — 1 = l — hu, thus at least for any 7 > 1/2 since > 1/2. 
Here and in what follows variables with the D index refer to the random bond case. 
Equation (|87D, compared with equation ( |69D leads to 

(l + /f^)(7-r^ + l) + 2< 1 + 27 ,1/2<7<1. (88) 
In the case of self-affine behaviour, using eq. (pOD, equation (^) gives: 

Hd>1 , l/2< 7 <1. (89) 

As before, for the 7 = 1 case the above inequalities are useless. In the case of 
self similar behaviour, inequalities for the fractal dimension di^D can analogously be 
deduced: 



din < 1 



, l/2< 7 <1 

H , equations 



and 



(90) 
OD together 



Because, in general < H < I and I < di < 1 
yield for 7 e (1/2,1): 

= 1, di^D = 1, i.e. the mean field values. 
The 7 = 1 case has been exactly solved |^ , and it can be mapped into a problem of 
a directed polymer in a random medium or equivalently a domain wall in a disordered 
two dimensional ferromagnet P2|, |53|, I 



for which an exact solution is known. 



The corresponding values for the exponents are: 



r = 7/5 , ip = 5/3 , H = 2/3 , di = l ,h = 3/5. (91) 

Disorder can be introduced in the system in another way, replacing the constant 
injection in each site of the lattice with a random, quenched, local injection, i.e. a 
spatial inhomogeneity in the rainfall in eq. (|l|). 

In order to do that, one can associate with each site i of the lattice a random variable 
Tj. The variables are chosen to be independent of each other, identically distributed 
and with mean (r^) = 1 . 

The accumulated areas must then satisfy, as in 
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A = 5Z '^id^j + ' (92) 
j 

in such a way that 

{1 if i is connected to j 
through upstream drainage directions /„„x 
orifj = i, ^^^^ 
otherwise . 

The minimum of the energy averaged over the "random-rainfall" will be denoted by 
Erri'y), and for a given value of 7 is given by 

E,,(7) = (min^A,({r,},T)^), (94) 

i 

where S denotes the set of all spanning trees T and {rj} denotes the whole set of 
random variables. As in (|85|) T* will be one of the trees for which the minimum of 
the energy is realized in the absence of randomness in the rainfall and for a given 
value of 7. Then 

Erril) = {^^T.M{r^hTr)<{T.M{r^},T*{^)y) 

i i 

i j i j 

i 

Thus 



Erril) < E{j) ~ min(L2, L^+^t) _ (gg) 

In this case, it is possible to bound the energy also from below with an argument 
analogous to that used in Section 4 for the homogeneous case. The detailed calculation 
is given in the appendix. Thus one can conclude that 

Errh) mm{L\ L^+^^) , (97) 
and all the results of Section 4 hold. 



6 Numerical Results: global minimum 

We have carried out extensive numerical investigations of OCN along two avenues: 
i) the search for the global minimum with a Metropolis algorithm for 7 = 1/2. 
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ii) the statistics over local minima for 7 = 1/2; strikingly these yield consistent 
but different values for the scaling exponents. By a local minimum, we mean a 
configuration (a spanning tree) of the network such that no link can be changed 
without increasing the energy. The global minimum is of course also a local minimum; 
but in the two cases we found different statistics, that is suggestive of a very rich 
structure of the energy-landscape. We will postpone the results concerning the latter 
subject to the next Section, focusing our attention only on the scaling properties of 
the global minimum. 

In the computer simulations we considered a square lattice with all sites on one 
side allowed to be an outlet fot the network. Periodic boundary conditions were 
chosen in the other direction. 

Once the simulation has been performed over the whole lattice, the basin with the 
biggest drained area is selected and only sites contained therein are used to calculate 
statistical quantities. Multiple outlets are allowed in order to minimize finite-size 
effects. The rainfall is assumed to be uniform over the whole lattice. The optimiza- 
tion method used is simulated annealing, in which a parameter T analogous to the 
temperature is introduced and lowered during the simulation. For each T value the 
system is relaxed in the following way: a new allowed configuration "near" the initial 
one ( "near" means one differing from the previous one only in one link) is randomly 
chosen; the dissipation energy of the new configuration is calculated and compared 
with the value of the old one. The new configuration is accepted with probability 1 
if AE is negative, and with probability exp[— ^] otherwise. 
In short a sketch of the algorithm is: 

i) Generation of a random initial configuration: 

Starting from a given, fixed tree we obtain a random initial configuration 
running steps ii) and in) several times. In this way only changes preserving 
the spanning loopless structure are accepted. 

ii) Random changes of the configuration: 

We select randomly one site i and then again randomly one of the "free" 
(not yet belonging to the tree) links connected with i, if any. 

Hi) Geometrical controls: 

The absence of a loop in the new configuration is checked. If a loop is 
present, step ii) is repeated. 

iv) Energetic control: 

The change, AE, in the dissipation energy dissipation is calculated. If it 
is negative we go on to step v). Otherwise, a random number p uniformly 
distributed in the interval [0, 1] is generated and compared with exp[— ^]: 
if exp[— ^] < p we go on to step v) , otherwise the change is rejected and 
we go back to step ii). 
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v) Recalculation of changed quantities: 

All variables involved in the change are updated. 

vi) Lowering of the T parameter: 

In each cycle, the T parameter is decreased with the following rule: at 
the n-th cycle T is given by T{n) = a"'T{0), where a is a parameter very 
close to 1 (we choose a = 0.986) and T(0) is a suitable chosen constant. 

After step i), steps u)-v) are repeated many times, say A^. Then we go to step vi) 
in which the T parameter is lowered and the entire algorithm (except step i) ) is 
repeated. The annealing process stops when T reaches very low values (~ 10~^). 

The simulations have been repeated varying the initial configuration for a size 
L = 128. The statistical quantities do not depend on the initial data. The integrated 
probability distributions for the accumulated areas and mainstream lengths averaged 
over 10 trials are shown in Fig. 5 and Fig.6 and give r = 1.50±0.02 and ip = 2.00±0.02 
where the error is estimated as the root mean square root over the ten trials. 
These results are in perfect agreement with equation (p9D and confirm that the ana- 
lytical results hold when 7 = 1/2. 



7 Numerical Results: local minima 



Homogeneous OCN yield results in good agreement with experimental data on rivers 
when statistics based on local minima is calculated. This suggested the concept 
of feasible optimality according to which nature is "unable" to reach the true ground 
state when complex systems are involved in optimization problems. The optimization 
just stops when one of the local minima is reached. Within this framework, the 
scaling properties of real rivers should be reproduced considering the ensemble of 
local minima. 

We have performed the search for local minima with an algorithm equivalent 
to a T = Metropolis scheme, i.e. one in which new configurations are accepted 
only if energetically favourable. The simulation has been repeated 40 times, starting 
with different, randomly chosen, initial data and varying the size L of the system: 
L = 32,64,128 . 

The values obtained for the characteristic exponents of the probability distribution 
considered above, r and ip are in very good agreement with experimental data. The 
distributions obtained starting from different initial condition do not substantially 
differ from one another, all yielding the same value for the exponents. 

The statistical quantities thus seem to be robust with respect to variations of 
initial conditions, showing the self-averaging of the scaling parameters. 

Results obtained averaging over 40 local minima are shown in Figs. 7 and 8 and 
give r = 1.45 ± 0.02 and = 1.82 ± 0.02. The results are consistent with the scaling 
relations. A collapse test has been done to verify the consistency of numerical values 
of the exponents with the finite size scaling hypothesis (Fig.9). 



21 



8 Summary and Conclusions 

In this paper, we liave studied the Optimal Channel Network for the drainage basin 
of a river. 

Within the framework of the finite size scaling hypothesis for the distributions of 
accumulated areas and mainstream lengths, we deduced the exact scaling behaviour 
for the tree (drainage network) for which the absolute minimum of dissipated energy 
is realized in both the homogeneous case and in the presence of randomness. The 
scaling exponents in the homogeneous case are found to be the mean field ones and 
differ from those measured in real rivers. 

Numerical results were obtained both for the statistics of the global minimum (con- 
firming the analytical predictions) and for local minima. The statistical exponents 
characterizing the local minima definitely differ from the mean field ones. They seems 
to be in a new distinct universality class and in agreement with data from real rivers. 
This suggests that real rivers, during their evolution, do not visit all of configuration 
space but are soon trapped in a metastable state, i.e. a local minima of the dissipated 
energy. 

The authors wish to thank A. Rinaldo and M. Cieplak for many enlightening discus- 
sions. This work was supported by NASA, NATO, NSF, The Petroleum Research 
Fund administered by the American Chemical Society and The Center for Academic 
Computing at Penn State. 



9 Appendix 

Using the notation of Section 4 we denote by P„ the lines orthogonal to the diagonal 
passing through the outlet, enumerated starting from the corner farthest from the 
outlet with 0. It will be useful to associate each line with n < L — 1 with a 

Vn = V^2L-l-n) • (98) 

In what follows we will choose the {rj} to be independent random variables with 
values in the interval [0, 2] and such that (rj) = 1. We then proceed to the first step 
of the proof as in Section 4, observing that the sum over all the sites in equation ( ^4]) 
defining the energy can be performed in two steps: 



Erril) > {mm J2A{{r,},Tr) 

i 

L—1 / n 2L—l—n \ -y 

> E((EE-^+ E Eri)), (99) 

n=0 ^k=l i 6Dfc k=l i ^ 
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where in the last step the equahty holds only for directed trees. 

The last expression can be written in a more convenient form introducing 



otherwise . 

Then 



^12 i eVk, with k <n 



Err > ^^2^L^^f^LJl^ifEJ^Y) 



n=0 

L-1 



2L2 



> 27r27 ^^=i L^i ev^ H'ty ii _ 27 T- 1+27 



(101) 



where the last inequality follows on observing that 



2L-l-n. 



k=l i eCfc 



Thus, for Ti < 2, the argument between the brackets is less than or equal to one. 

Furthermore, for any x and 7 G [0, 1], > x. 

Equation (|101|) together with equation (|96|) of Section 5 gives 



27^1+27 < ^^^(^) < ^(^) _ ^1- 



-27 



(102) 



and thus 



(103) 
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{a)9 ~ _ 

As noted in [^ the mean accumulated area for a border site in a basin is just 
the number of sites in the bulk divided by the number of sites in the border. 
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Figure 1: (l.a) The basin of Fella River in the north-east of Italy reconstructed from 
a Digital Elevation Map. (l.b) A lattice river basin of size L = 5. In each site i 
the value of the cumulated area Ai is displayed. The darkest line represent the main 
stream of the entire basin. 

Figure 2: (2. a) The Peano basin at iteration step T = 0, T = l, T = 2 and T = 3, 
with the accumulated areas displayed. (2.b) To obtain Peano basin at time step T+1, 
one keeps the basin at time step T, cuts the outlet and joins four copies of what is 
obtained as illustrated in the figure. In this way, the recursion relation can be easily 
understood. 

Figure 3: The renormalization group argument for the Peano basin. B-sites die under 
decimation. 

Figure 4: Each dashed line divides the lattice in two parts. J2iev is at least equal 
at the number of sites contained in the part of the lattice with border T> and without 
the outlet. 

Figure 5: The distribution P{A > a) versus a for a basin of linear size L = 128. Such 
a distribution has been obtained by mean of a Metropolis-like algorithm (attempting 
to seek the global minimum) and averaged over ten samples. The slope displayed is 
1 - T = -0.50. 

Figure 6: The distribution n(/ > A) versus Iq for the same conditions as in Fig. 5. 
The slope 1 --0 = -1.00 

Figure 7: P{A > a) versus a for three different sizes of the basin: L = 32, L = 64 and 
L = 128. Local minima are considered here, the distribution is obteined averaging 
over 40 samples. The slope is 1 — r = —0.45 

Figure 8: The same as in Fig. 7 for the upstream length distribution. The slope is 
1 - -0 = -0.82 

Figure 9: The result of the collapse test for the accumulated area distribution in the 
case of the local minima. 
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